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PERTURBATIONS OF 
NON-RESONANT SATELLITE 

ORBITS DUE TO A ROTATING EARTH 
by 

Alan Mueller 


1.0 INTRODUCTION 


The dominant perturbations of the motion of a satellite 
near the earth are due to the non-symmetrical gravitational 
field and the atmospheric drag. The gravitational field may 
be divided in two classes: terms independent of time (zonal 

harmonics) and terms which depend explicitly on time. A com- 
plete first order solution for satellite motion perturbed by 
the second zonal harmonic has been developed by Scheifele (ref- 
erence 1). This solution has been rewritten in the non-singu- 
lar PS(}> elements by Bond (reference 2). In references 3 
and 4, the perturbations due to drag are developed in PS4> 
elements and added to this theory. In reference 5, the 

long period perturbations due to the additional zonal harmon- 
ics are included in the theory. The tesseral harmonics have 
yet to be treated under this unified theory and are the topic 
of this report. 

The perturbations due to the tesseral harmonics can be 
placed in four categories: 

1) Short period perturbations with a magnitude of 


t 


A first order solution here implies that the solution has 
a periodic, error of second order and a secular error of 
third order. 



“ 6 “ 
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about J| 

2) Intermediate period perturbations with a magnitude 

of between Jj Jj . 

3) Along track secular perturbations induced by the 
periodic perturbations in the mean motion. 

4) Resonant perturbations. 

For near earth satellites, the atmospheric drag perturba- 
tion continually pulls the satellite in and out of the differ- 
ent long period resonant frequencies. The result is that the 
resonances never become apparent and may be neglected. 

Since the theory has been developed only to first 

order, the second order short period tesseral perturbations 
may be neglected. For the same rer.son, the intermediate pe- 
riod perturbations should be included if first order accuracy 
is to be maintained. The tesseral hamwnics have no true sec- 
ular perturbation but the periodicities in the mean motion 
induce a secular perturbation in the mean anomaly. This sec- 
ular perturbation may be determined by simply using the aver- 
age mean motion instead of the osculating mean motion. Graf 
(reference 6) finds the average mean motion in a numerical 
manner. The numerical studies in that reference verify the 
assumption that use of the average mean motion does account 
for the apparent secular trend in the mean anomaly. Although 
the results are good, use of a numerical method would be in- 
consistent with the idea of having a completely analyt cal 
theory. 

To complete the solution of the motion of a near earth 
satellite, the averaged mean motion and the intermediate pe- 
riod perturbations need to be found in a completely analytical 
manner . 

Since the previous PS0 theory has applied Von Zeipel's 
solution technique, it seems natural to return to this method 
for a solution of the tesseral perturbation. As in the previous 


theory, the solution will he found first in the singular DSiff 
eleoaents and then rewritten in the PS^ elements to remove 
singularities. The notation used in the development is de- 
scribed in Appendix A. 
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2.0 FORMULATION OF THE PROBLEM 

Only the gravitational field will be considered since 
its Interaction with the drag effects is small. The DS0 
hamiltonian for the gravitational potential is assumed to 
have the form 


where 


F = Fq + eFj + e^F^ 


( 1 ) 


F a * 

0 M 


( two-body ) 


( 2 ) 




(Jo) 


(3) 


F = F + F 
^2 Z T 


(Zonals and Tesserals) (4) 


Vi 


R*' 


' p° (^) 


' q[-‘ fr'} r"-' ■" 


(5) 


U 


n=2 m=l r 


C cos m( X-6 ) + 
n , m 


+ S sin m(A-G) 
n , m 


( 6 ) 


3 

f = - J. u R' 
2 2 


(7) 
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P® are the associated Legendre polynomials; Is the mean 

equatorial radius; J_ , C , S are the geopotential 

^ 4 n»n n,iD 

coefficients; X is the longitude of the satellite and 9 is 
the Greenwich hour angle. 

The hamiltonian may be divided in zonal and tesseral 

harmonics. The functional dependence of and F^ is 

Fz = Fz(4»,g.B) . (8) 

F.J, = F^((|>.g,h,a)^il'’’,6) . (9) 


The canonical time element I will always appear, in the 
problem, premultiplied by the rotation rate of the earth 

0 
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3.0 VON ZEIPEL'S SOLUTION TECHNIQUE 

As was stated earlier, the Von Zeipel solution method is 
to be used to find the average mean motion and to eliminate 
the periodicities due to the tesseral harmonics. The genera- 
ting function S is to be used for the elimination of the 
short period function of and the intermediate period func- 
tion of . The equations governing the transformation in 

the DS4) space are 

as as 

o' = — , e = — . (10) 

ae' 3 a 


As in the hamiltonian, S is assumed to be of the form 


S » Sg t-ZcS, + 


( 11 ) 


s, - ^ (12) 

where S.^. is a periodic function of (}> and and 

and are periodic’ functions of (J> only. The development 

of Sj can be found In reference 1 and therefore the discus- 
sion will be restricted to the development of $2 

The theory of reference.s 1 and 2 neglected S 2 to main- 
tain a first order accuracy. But this assumes that S 2 is 

order 0(1). This is true for zonal perturbat icjtis but the 
tesseral perturbations result in an S 2 which may become 

larger than 0(1). This will be demonstrated in the following 
analysis. The necessary Von Zeipel equations found from ref- 
erence 1 are 


F, 


0 


: (4>* .1/ ) = Fq (iJ)' ,1/ ) 


(13) 
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dF as, 

- p (3'.0.g) + f; (B’.r) (14) 

3$ 3^ 


: !!o !fi ^ !!o 

- --T' 

^ a Fq 3Sj aSj 

34* 34 3L* 3£ 

2 ^ ^ 

^ 38^3Bj aa^ 

A i 

kM 3B’ 3a^ 6 

'* 

3P; as, 

► L L 

:f 3a^ 33- 

-F^CS'.^.g)- (15) 

- F^(6’,4,g.h,a)^JO 

+ FJ(B’ ,g) 

• 

Assuming S 2 has the form given by 12 

, then the expres- 

sions required to find 

and are given by 

^''0 *®2 . ‘ V V 

3^Pq 3S, 

3S, . 

^ 4HH# 

34>'34 2|(=i 

3a^ 

3a^ 

k»l 

.y 

k“ 1 k 

^^1 - 

( 16) 

ssic 


- F^{B* .(j ,g) - F’(6' .g) 






3L’ ?a 


( 17) 


If is to be used in eliminating th€‘ intermediate period 

terms then it is necessary to keep only those terms of order 
greater than 0(1) . If is to be used in finding the 

f 

♦4 
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mean energy, only those terms which are functions of t are 
considered. With these qualifications a much simpler set of 
equations can be used to find 8^ i 

Sj - 0 , (18) 


3Fo as^ 3 Fq as^ 

+ - F (B'.<^,g,h,w 1) 

a^' H aL* at ^ ® 


(19) 


S 2 is set to zero because the right hand sides of equation 

(16) are of order 0(1) and are Independent of I 

Thus the transformation due to the short and intermediate 
terms is given by 


where 


3S, 

0 + e + £ 

2 I 

as' 

38’ 

as 

6’ + £ i- + 

as_ 

e* * 

3a 

da 

0(1) 


element L 

is the total 


( 20 ) 


( 21 ) 


lated to the mean motion. The mean total energy L' can be 
then related to the instantaneous total energy by 


as 

L - L’ 4 . C* — - 
di 


( 22 ) 
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where 


St i 


0 ( 1 ) 


Remarks: 

If the Von Zelpel's solution method had been applied to 
the classical Delaunay variables a much more complicated set 
of equations would appear. The average mean motion could also 
be found by differentiating the classical determining function 
S with respect to the mean anomaly. But since (or F^) 

Is a function of the mean anomaly, S 2 could not be considered 

zero since it would also be a function of the mean anomaly. 

Also the expression for the mean energy would be more complex 
and difficult to solve since terms of 0(e) would appear. All 
the usual ccmipllcations cf the second order theory would appear 
just to maintain a first order accuracy. But with the DS4> 
elements the equations for finding the generating function re- 
main concise and relatively easy to solve. 

The mean energy is found through a determining function 
which is defined by a partial differential equation with the 
anomaly of the satellite 0 and the changing hour angle 
u. t . This expression is very similar in form to the 

w 

part al differential equation found by Graf (releronco 6) by 
applying the "Method of Averages" to the KS differential 
equations. Maybe it is not so peculiar that the Method of 
Averages and Von Zeipel’s equations result in similar formulas. 

If can be written in a fourler series 01 the form 

T 




k J 


Cj^^(S'.h.g) .f08 (k^ - jujg,i) 4- 


(B'.h.g) ain (k(J» - 


( 23) 


then the solution for in equation (19) can be easily shown 
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to be 




( 24 ) 


wh^re V is the ratio of the frequency of the earth's rotation 
to the revolution frequency of the satellite. 


U), 


O 


3L 


= w 


®(2L)^" 


(25) 


For near earth satellites this ratio is about one to sixteen, 


1 - 

V % — 
16 


(26) 


Therefore when k = 0 and for 



1 

the factor 

k-v.1 


becomes on the order of O(^) • These terms result in the 

intermediate period perturbations which are the order of 

o(4> • 


In order to complete the solution, the tesseral Hamilton- 
ian F must be expressed as a fourier series of the form 
T 

given by equation (23). As before the solution will then bo 
expressea in the non-singular PS(|i elements. 
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4.0 EXPANSION OF THE TIME DEPENDENT GEOPOTENTIAL 


The definition of as given by equation (6) is 


n 


F, - - E r r V 

T 2 " n,m 

n=»2 m=l 


(27) 


where 


r" ^ r 

Jy = p“ (|) C r-os m (X-e) + 

n,m y."~^ " ^ |^n,ra 


+ S sin m (X-9) 
n , m 


(28) 


We desire an expansion of the above expression in a 
fourier series of the canonical DS<|) angular variables (J> , 
g, h, and I . The expansion will, in many respects, be sim- 
ilar to Kaula's expansion of the disturbing function in the 
classical angular variables (reference 7). The important dif- 
ference is that expansions will now be carried out in the true 
anomaly instead of the mean anomaly. 

As in Kaula, the potential may be expanded using the in- 

clination F ' to a form 
ntnp 



( 29) 


To maintain the notation of previous authors, this author 
has decided to keep the notation p as an index in the 
inclination function. For distinction, the semi-latus rec- 
tum will now appear as an italic "p" 






irtiere 


4>nnp * (n-2p) (<|>+g) + in (h-0) 


(30) 


nm 


nm 


-S 


nm 


n-m even 


n-m odd 


(31) 


B 

nm 


nm 


nm 


n-m even 


n-m odd 


(32) 


The recursive relations for the inclination function have been 
derived by Giacaglia (reference 8). These relations and others 
appear in Appendix B. 

Another transformation is necessary to replace the powers 
of r as a fourier series in 0 : 



aos 

sin 


nmp 



E 


j*l-n 



cos 

sin 


\jj 

nmpj 


(33 ) 


"^nrapj “ (n-2p+j) (p + (n-2p) g + m (h-6) 


(34) 


Note that this Is a finite expansion with an eccentricity 
function (e) . The development of g"" is given in 

Appendix C. 

A final transformation is required to place the disturb- 
ing function in the form of equation (23). In the expression 
of the angle ip , the hour angle may be written as 


e = 0^ + o),-t 

o © 


( 35 ) 
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where 6^ is the hour angle at time t « 0 . Time is now 

a dependent variable^^ l^ t^% DS^ system given by 


t » «. + 


(2L) Va 


[m (♦) - *] 


( 36 ) 


where 

M(4») = E - e sin E . (37) 

Replacing equations (35) and (36) into equation (34) one finds 


’^nmpj " (n-2p+j) 4) + mv (c{i-M) + (n-2p) g + 


+ m ( ) 


Thus we define a new function N (v,e) which gives 

Iv 


( 38) 


aos 

^nmpj 

k 

cos 

sin 

k 


sin 




nmp j k 


(39) 


where 


~ (n-2p+j+k) ({) + (n-2p) g + m (h-(;j Jl-e ) . (40) 

nmp j K © o 

The development of N™ can be found in Appendix D, 

Accumulating the results we arrive at the final expres- 
sion for the disturbing function ; 
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tt n n* 1 * 

n«2 m*l p«o j“l-n k“-« 


(41) 


where one element is given by 


B sin i|; 

nmpjK nmpjk | nm nmpjk nm nropjk 


(42) 


and 


R 


n 

® F n“ 


nmpjk ^ ^ 


(43) 


Thus one element of the determining function defined by 

equation (19) becomes 


-J 


Tntnp j k 


nmp.i k 
(n-2p+j+k-vm) 


A sin 4> 


nm 


nmpjk 


(44) 


- B SOS i|j 

nm nmpjk 


The mean energy can then be found from equations (44) and (22); 


q nmpjk nmpjk 


where AL ,, is expressed as 
nmpjk ^ 
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nmp jk 


nuD. 


iL?.PiK 


(n-2p+j+k 


-vm) I 


A ao8 + 

nm nmpjk 


( 46 ) 


+ B sin ip 

nm nmpjk 


The elimination of the intermediate period terms may also be 
found by determining the partial derivatives; of each element 
of S„ in which 

Tnmp j k 


k = 2p - n - j 


( 47 ) 
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5.0 EXPRESSING THE GEOPOTENTIAL IN PS^ ELEMMTS « 


In the expansion of the geopotential, the angle 4'nnpjk 

becomes undefined for small eccentricities and inclinations. 
The fourier expansion of equations (41), (42) and (43) should 
be transformed into a fourier series in the true longitude Oj 

and the angle w-Jl+O , and polynomials of the well defined 

w O 

functions 


ri * e ein (g+h) 
C = e 008 (g+h) 


6 “ s * 8in h 

g = S • 003 h 


(48) 


With the definition of the true longitude » 0 + g + h 
the angle ’J'n^pjk 


= 9 “ (j+k)(g+h) + (m-n+2p) h , (49) 

nmpjk nmpjk 


0 = (n-2p+j+k) a, - m(a)^S,+6 ) . (50) 

nmpjk u ^ I 0 o 

Keeping the notation of Giacaglia (reference 9) we make new 
definitions 


q = j + k 


a = m - n + 2p 


so that tj; becomes 

nmp i k 


\p =0 “ q(K+h) + ah 

nmpjk nmpjk 


(51) 

(52) 


(53) 


To maintain the notation of previous authors, this author 
has decided to keep the notation of q as an index in the 
expansion. Please observe the definition of q given in 
Appendix A. 


t 
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By defining 


^ (g+h) 6 q, •* 8^“^ ein a h , 

Cq “ <j(?e q (g+h) aoa o. In 


R = c 6 + n 6 

aq q a q a 


I = C 5 - n 3 

aq q a q a 


the expressions for the functions of \b become 

^ ’^nmpjk 


el^igl^l ^lJ = R CC7S .6 .. + 

nmpjk aq nmpjk 


+ I stn 0 

aq nmpjk 


•u'^R^ sin Q ., + 
nmpjk aq nmpjk 


- I aos 0 

aq nmpjk 


Recursive relations exist for the definitions of 
tion (54). They are 


n = n , i; , + c , n , 

q q - 1 1 q- 1 1 


S " ^q-1 "l 


^ ^ 

-q q 


^.q = S 


= Q P 2 


Hi = - Q a. ^ , 


(54) 

(55) 

(56) 

(57) 

(58) 

equa- 

(59) 

(60) 
(61) 


(62) 
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«»■ 

«a-: 

Si * Vi «i 

(63) 

■ 


«1 - • 

(64) 

6 S3 

-a 

-6 

a 

■ ®-a " 

(65) 

“ 

ppj 


(66) 


With the expressions (57) and (58), one is now able to 
write each element of the disturbing function given by equation 
(42) entirely in non-singular coordinates: 


= J 


nm 




A R - B I ) e . 

nmpjk 


pjk nmpjk | nm aq nm aq 
+ (A I +B R ) St n 6 


( 67 ) 


nm aq nm aq 


nmpjk 


J 

nmpj 


k 



e 


n - 1 

P 


F 

nmp 


gn-l 


g( i j l + lit|-| j-'kl ) 


( 68 ) 


F 

nmp 



( 69 ) 



( 70 ) 




( 71 ) 


The expressions of the barred values are found from the 
relations of the unbarred values. These barred expressions 
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are free of singularities. The relations for f 
and may be found in Appendices B, C and 

respectively. 


nmp * 
D. 
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6.0 CONCLUSIONS 


Orbit perturbations due to the time dependent harmonics 
of the gravitational field of the earth have been studied for 
the case of near-earth satellites. Since for these orbits, the 
atmospheric drag perturbations never permits the satellite to 
remain long in a resonant frequency, the resonant perturbations 
have been neglected. Because the short period perturbations 
due to the tesseral harmonics are of second order, these terms 
have also been neglected . 

The solution for the satellite motion which includes the 
intermediate period perturbations of between first and second 
order have been found by using the Von Zeipel solution method 
with a general, recursive expansion of the geopotential in the 
non-singular PS<J) elements. Also, the along track secular 
perturbations induced by the periodic perturbations in the 
mean motion are eliminated by the computation of the average 
energy using the same theory. 

This theory has been implemented in an operational can- 
puter program ASOP (reference 12) and numerical experiements 
verify the expected accuracy. 
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APPEHDIX A 
NOTATION 


The definitions of the DS^ elements are as follows : 
The angular variables 

ttj ■ 0 true anomaly 

a- ■ g argument of pericontor 

(Al) 

= h ascending node 

a, * Z time element 

k 


The action variables 

fij “ ^ related to two-body energy 

B. > G angular momentum magnitude 

(A2) 

= H Z component of angular momentum 
6^ = L the total energy 

These may be transformed to the canonical PS(j) elements: 

Oj=0 + g + h Pj=<J> 

O 2 = - ^2( $-G) ain (g+h) P 2 ” \/2( <I>-G ) i‘o& (g+h) 

^ (A3) 

= - >/2(G-H)' Bin h = 72(0-17) h 


Abbreviations used in the text 



are 

(semi-latUK rt’ctum) 
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H 

C ■ COB I ■ — 

G 


(eccentricity) 


(1 is the inclination) 


Q 


P = 




2p 


+ G - 


y2(G+H)' 


2G 



APPENDIX B - INCLINATION FUNCTION 
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APPENDIX B 
INCLINATION FUNCTION 


Recursive relations for the Inclination Function F 

nmp 

have been derived by Giacaglia (reference 8), The first re- 
currence relation given in that reference is simple and valid 
for all n , m , p except when n = m : 


F 

nmp 


^ , (2n-l)s r I 

jn-m+l j, -F - 

2(n-m) L n-l.”.P-lJ 



2 , m, p 


(Bl) 


The other relations derived in that reference were not 
suitable because they had singularities for zero inclinations. 
A new recurrence relation derived through induction by this 
author completes the definition of the inclination function: 


2n-l 


“1 


F = <sMl+c) F , , , + 

nnp 2 I n-l,n-l,p-l 


+ (1+c) F 


n- 1 , n-1 , p I 


(B2) 


In all the recurrence relations 


F = 0 
nmp 


for p < 0 and p > n 
m < 0 m > n 


(B3) 
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functlon is 


a ■ in - ,n + 2p . (B4) 


From equation (Bl) one finds the relation for F is 

nnp 


F 

nmp 


(- 1 ) 


n-m+ 1 


2n-l 


2(n-m) 


a F 


n- 1 ,m, p 


- b F 


n-1 ,m,p-l 



(n-l+m) 

(n-m) 


F 


n-2 ,m , p- 1 


(B5) 


The values of a and b are found from the conditions 


n >m+2p 

n<m+2p 

n=m+2p 

a = 1 

s^ 

s ^ 

b = s ^ 

1 

s ^ 


From equation (B2) the definition of F becomes 
^ nnp 


2n-l 


nnp 


(1+c)'^ , , , + 

I n~l,n-l,p-l 


+ (1+c) F 


n-l ,n- 1 , p 


(B6) 
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As before the following applies 

? = 0 for p < 0 and p > n 

niop 

m > 0 and m > n 

The starting value of the recurrence relation is 


DOT 


(B7) 


1 


(B8) 
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APPENDIX C 
ECCENTRICITY FUNCTION 


The eccentricity function 
relation 



is defined by the 



P 

where — = 1 + e cos t}i 

r 


The binomial expansion of 




becomes 


cos (s-2t) 


(Cl) 


(C2) 


Using equation (C2) and the definition of equation (Cl) 
the expression for g" becomes 



for j ^ 0 


(C3) 


-44- 


Recurrence relations exist for the binomial coefficient 


\k) (n-k)!k! y k j 


(C4) 


The expression for g" becomes 




-U = 


-j 


0 - 


=E 

kHJl 


.k-j 


(v)(j 


(C5) 
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APPENDIX D 
THE TIME EXPANSION 


It Is desired to make the following expansion 


mv (4>-M) = 


■ ■ 

‘L.K . • 

k Btn 

L J 


To develop this expansion, a complex notation is adopted 
which has the same meaning 


exp [imv (4)-M)] = N® exp (ik(j)) 


where is defined by 


1 

= — I oxp [imv (4>-M) - ik({)l dcj) 

2-n J '■ 


As In the development of Hansen coefficients (re*ference 10) 
some abbreviations are useful 


y/ 1 +e' 1+6 


z = exp iE 


(D5) 


The relatioa between the true and eccentric anomaly 



tan - ^ ■ tAn ^ E 

2 /W ‘ 2 


(D6) 


can be written in another form by equations (D4) and D5) 
as 

exp i(J) = = z (1- f) (l-3z) . (07) 

i-3z ^ 


Similarly the relation between mean and eccentric 
anomalies 


M = E - e sin E 


(D8) 


may be written as 


exp iM = z • 



(D9) 


By differentiating equation (D7) , the expression for d4) 
becomes 


( 1 - 6 ^) 


(1-Bz)"^ z'^ 


( 1 - 


f>- 


dz 


d4> = 


i 


(DIO) 
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Tbus the integrand in equation (D3) becomes 


exp i 


1 \ z 



9 \"(k+l“»v) 


. exp[^ ( 4 )] 


dz 


.-1-k 


(Dll) 


Tne powers of (1-Sz) and 1-^ can be expanded with 

z 

small parameter 6 


f = (1-gz) 
s 


k-mV- 1 


df 


- = (k-mv-1) 


dz 


d'f 


dz' 


® (k-mv~l )( k+mv-2 ) 


(D12) 


s 


s - 1 


n [k-(mv+j+l)] 
j*0 


d^f (z=0) 
s 

J s 

dz 


s- 1 


= -(B)*' 77 [k-(mv+j+l)j 


j*0 
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By Taylor's expansion 


f - 1 + 

8 


4*f(z-0) 1 


B¥t dz si 


or 


E s 


g«0 


where 


(- 1 ) 


8 8-1 


a = 1 a 

° ® s! 


77 [k-(mv+j+l)] 


J-0 


or 


a 


( k-itiv-s \ 


Similarly 


= (l-6w)’ where w “ - 

z 


df 


(k-mv+l) 


dw 


d^f 


dw' 


t = B^(l-Bw) 


(D13) 


- (k-mV+3) 


( k-mv+l )(k-mv+2 ) 


(D14) 
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t c-l 

tJs m (k-mv+l+j) 


d^f^(w»0) 


dw‘ 


t-l 

n <k-mv+l+J) 
j-0 


Therefore: 


f 


t 





or 


go 





(D15) 


wher<* 


b 


o 


= 1 


b = — 77 (k-mv+l+j) 

t j = 0 


or 


( 


k-mv+t 



With equations (D13) and (D15) the integral becomes 
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But the Bessel function is defined as 


J (a) 
n 


1 

2iri 




so that the final form becomes 


(D17) 


N” -•* (1-BM 




8 + t 


k-8 + t 


(mve) 


s«0 t»0 


(D18) 


where 




( k-mv~s ) 


a = - 
s 


8-1 


(D19) 


b = 1 

o 


( k-mv+t ) 

b b , 

t . t- 1 


(D20) 


The barred function N™ can be found from its definition 

k 

and equation (D18) 


m 

N™ = 

e ' 


(D21) 


If one notes the definition of 0 in equation (D4) , 
one finds 

00 OD 


K - (1-8’ EE- 


CmjJLdt yj. /TTT\8 + t 


20- ,|k-9 + t(t , . fnnn\ 

e (mv) *^k-s + t (D22) 
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The value of j is determined by the following conditions 



(D23) 


Note that o can never be negative. 

The function J (a) is defined from the Bessel function 
n 

as 


Jn<a) 


J^(a) 



(D24) 


It can be found from following recursive relation 


n- 1 


= 2n J 


- a' 


n + 1 


( D25 ) 


with starting values 



(D26) 
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The expressions tor N™ were determined in another 

manner by Bond (reference 11) up to order e^ . These ex- 
pressions are given here in the barred form. 



(D27) 



